The lnRR_func function is here used to calculate a log response ratio (lnRR) adjusted for small sample sizes. In addition, this formula accounts for correlated samples. For more details, see Doncaster and Spake (2018) Correction for bias in meta-analysis of little-replicated studies. Methods in Ecology and Evolution; 9:634-644
# packages
library(tidyverse)
library(googlesheets4)
library(here)
library(metafor)
library(metaAidR)
library(orchaRd)
library(ape)
library(clubSandwich)
library(metaAidR)
library(patchwork)
library(emmeans)
library(kableExtra)
library(GGally)
# custom functions - here
# this is to get small sample size corrected lnRR
# also, the correlated samples - MS should put formulas of what we used and assumptions
# averaged
lnRR_func <- function(Mc, Nc, Me, Ne, aCV2c, aCV2e, rho = 0.8){
lnRR <- log(Me/Mc) +
0.5 * ((aCV2e/Ne) - (aCV2c/Nc))
var_lnRR <- (aCV2c/Nc) + (aCV2e/Ne) +
rho*(sqrt(aCV2c)*sqrt(aCV2e)/(Nc+Ne))
data.frame(lnRR,var_lnRR)
}
# Mc: Concentration of PFAS of the raw (control) sample
# Nc: Sample size of the raw (control) sample
# Me: Concentration of PFAS of the cooked (experimental) sample
# Ne: Sample size of the cooked (experimental) sample
# aCV2c: Mean coefficient of variation of the raw (control) samples
# aCV2e: Mean coefficient of variation of the cooked (experimental) samples## # A tibble: 512 x 47
## Study_ID Author_year Publication_year Country_firstAu~ Effect_ID
## <chr> <chr> <dbl> <chr> <chr>
## 1 F001 Alves_2017 2017 Portugal E001
## 2 F001 Alves_2017 2017 Portugal E002
## 3 F002 Barbosa_20~ 2018 Portugal E003
## 4 F002 Barbosa_20~ 2018 Portugal E004
## 5 F002 Barbosa_20~ 2018 Portugal E005
## 6 F002 Barbosa_20~ 2018 Portugal E006
## 7 F002 Barbosa_20~ 2018 Portugal E007
## 8 F002 Barbosa_20~ 2018 Portugal E008
## 9 F002 Barbosa_20~ 2018 Portugal E009
## 10 F002 Barbosa_20~ 2018 Portugal E010
## # ... with 502 more rows, and 42 more variables: Species_common <chr>,
## # Species_Scientific <chr>, Invertebrate_vertebrate <chr>,
## # Moisture_loss_in_percent <dbl>, PFAS_type <chr>, PFAS_carbon_chain <dbl>,
## # linear_total <chr>, Choice_of_9 <chr>, Cooking_method <chr>,
## # Cooking_Category <chr>, Comments_cooking <chr>,
## # Temperature_in_Celsius <dbl>, Length_cooking_time_in_s <dbl>, Water <chr>,
## # Oil <chr>, Oil_type <chr>, Volume_liquid_ml <dbl>, Cohort_ID <chr>,
## # Cohort_comment <chr>, Nc <dbl>, Pooled_Nc <dbl>, Unit_PFAS_conc <chr>,
## # Mc <dbl>, Mc_comment <chr>, Sc <dbl>, sd <chr>,
## # Sc_technical_biological <chr>, Ne <dbl>, Pooled_Ne <dbl>, Me <dbl>,
## # Me_comment <chr>, Se <dbl>, Se_technical_biological <chr>,
## # If_technical_how_many <dbl>, Unit_LOD_LOQ <chr>, LOD <chr>, LOQ <chr>,
## # Design <chr>, DataSource <chr>, Raw_data_provided <chr>,
## # General_comments <chr>, checked <chr>
# creating SD for just biological
# not quite sure why if_else does not work but ifesle does (handling NA???)
dat <- processed_data %>% mutate(SDc = ifelse(Sc_technical_biological == "biological", Sc, NA), # Calculate the SD of biological replicates for control samples
SDe = ifelse(Se_technical_biological == "biological", Se, NA)) # Calculate the SD of biological replicates for experimental samplesThe phylogenetic tree was generated in the tree_cooked_fish_MA.Rmd document
tree <- read.tree(here("data", "plot_cooked_fish_MA.tre")) # Import phylogenetic tree (see tree_cooked_fish_MA.Rmd for more details)
tree <- compute.brlen(tree) # Generate branch lengths
cor_tree <- vcv(tree,corr = T) # Generate phylogenetic variance-covariance matrix
dat$Phylogeny <- str_replace(dat$Species_Scientific, " ", "_") # Add the `phylogeny` column to the data frame
colnames(cor_tree) %in% dat$Phylogeny # Check correspondence between tip names and data frame## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
The average coefficient of variation in PFAS concentration was calculated for each study and treatment, according to Doncaster and Spake (2018) Correction for bias in meta-analysis of little-replicated studies. Methods in Ecology and Evolution; 9:634-644. Then, these values were averaged across studies and used to calculate the lnRR corrected for small sample sizes (for formula, see the lnRR_func above)
aCV2 <- dat %>%
group_by(Study_ID) %>% # Group by study
summarise(CV2c = mean((SDc/Mc)^2, na.rm = T), # Calculate the squared coefficient of variation for control and experimental groups
CV2e = mean((SDe/Me)^2, na.rm = T)) %>%
ungroup() %>% # ungroup
summarise(aCV2c = mean(CV2c, na.rm = T), # Mean CV^2 for exp and control groups across studies
aCV2e = mean(CV2e, na.rm = T))
effect <- lnRR_func(Mc = dat$Mc,
Nc = dat$Nc,
Me = dat$Me,
Ne = dat$Ne,
aCV2c = aCV2[[1]],
aCV2e = aCV2[[2]],
rho = 0.8) # Calculate effect sizes
dat <- dat %>%
mutate(N_tilde = (Nc*Ne)/(Nc + Ne)) # Calculate the effective sample size
dat <- cbind(dat, effect) # Merge effect sizes with the data frame
VCV_lnRR <- make_VCV_matrix(dat, V = "var_lnRR", cluster = "Cohort_ID", obs = "Effect_ID", rho = 0.5) # Because some effect sizes share the same control, we generated a variance-covariance matrix to account for correlated errors (i.e. effectively dividing the weight of the correlated estimates by half)# mean
ggplot(dat, aes(x=lnRR))+ geom_histogram(fill = "salmon", col = "black", binwidth = 0.2) + theme_classic()# variance
ggplot(dat, aes(x=var_lnRR))+ geom_histogram(fill = "salmon", col = "black", binwidth = 0.05) + theme_classic()# log variance
ggplot(dat, aes(x=var_lnRR))+ geom_histogram(fill = "salmon", col = "black", binwidth = 0.05) + scale_x_log10()+theme_classic()dat %>%
summarise( # Calculate the number of effect sizes, studies and species for the main categorical variables
`Studies` = n_distinct(Study_ID),
`Species` = n_distinct(Species_common),
`PFAS type` = n_distinct(PFAS_type),
`Cohorts` = n_distinct(Cohort_ID),
`Effect sizes` = n_distinct(Effect_ID),
`Effect sizes (Oil-based)` = n_distinct(Effect_ID[Cooking_Category=="oil-based"]),
`Studies (Oil-based)` = n_distinct(Study_ID[Cooking_Category=="oil-based"]),
`Species (Oil-based)` = n_distinct(Species_common[Cooking_Category=="oil-based"]),
`Effect sizes (Water-based)` = n_distinct(Effect_ID[Cooking_Category=="water-based"]),
`Studies (Water-based)` = n_distinct(Study_ID[Cooking_Category=="water-based"]),
`Species (Water-based)` = n_distinct(Species_common[Cooking_Category=="water-based"]),
`Effect sizes (No liquid)` = n_distinct(Effect_ID[Cooking_Category=="No liquid"]),
`Studies (No liquid)` = n_distinct(Study_ID[Cooking_Category=="No liquid"]),
`Species (No liquid)` = n_distinct(Species_common[Cooking_Category=="No liquid"]),) -> table_sample_sizes
table_sample_sizes<-t(table_sample_sizes)
colnames(table_sample_sizes)<-"n (sample size)"
kable(table_sample_sizes) %>% kable_styling("striped", position="left")| n (sample size) | |
|---|---|
| Studies | 10 |
| Species | 39 |
| PFAS type | 18 |
| Cohorts | 153 |
| Effect sizes | 512 |
| Effect sizes (Oil-based) | 303 |
| Studies (Oil-based) | 7 |
| Species (Oil-based) | 28 |
| Effect sizes (Water-based) | 140 |
| Studies (Water-based) | 8 |
| Species (Water-based) | 23 |
| Effect sizes (No liquid) | 69 |
| Studies (No liquid) | 2 |
| Species (No liquid) | 14 |
kable(summary(dat), "html") %>% kable_styling("striped", position = "left") %>%
scroll_box(width = "100%", height = "500px")| Study_ID | Author_year | Publication_year | Country_firstAuthor | Effect_ID | Species_common | Species_Scientific | Invertebrate_vertebrate | Moisture_loss_in_percent | PFAS_type | PFAS_carbon_chain | linear_total | Choice_of_9 | Cooking_method | Cooking_Category | Comments_cooking | Temperature_in_Celsius | Length_cooking_time_in_s | Water | Oil | Oil_type | Volume_liquid_ml | Cohort_ID | Cohort_comment | Nc | Pooled_Nc | Unit_PFAS_conc | Mc | Mc_comment | Sc | sd | Sc_technical_biological | Ne | Pooled_Ne | Me | Me_comment | Se | Se_technical_biological | If_technical_how_many | Unit_LOD_LOQ | LOD | LOQ | Design | DataSource | Raw_data_provided | General_comments | checked | SDc | SDe | Phylogeny | N_tilde | lnRR | var_lnRR | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Length:512 | Length:512 | Min. :2008 | Length:512 | Length:512 | Length:512 | Length:512 | Length:512 | Min. : 6.77 | Length:512 | Min. : 3.000 | Length:512 | Length:512 | Length:512 | Length:512 | Length:512 | Min. : 75.0 | Min. : 120.0 | Length:512 | Length:512 | Length:512 | Min. : 5 | Length:512 | Length:512 | Min. : 1.00 | Min. :1.000 | Length:512 | Min. : 0.002 | Length:512 | Min. : 0.0010 | Length:512 | Length:512 | Min. : 1.00 | Min. :1.000 | Min. : 0.0020 | Length:512 | Min. : 0.000 | Length:512 | Min. :1.000 | Length:512 | Length:512 | Length:512 | Length:512 | Length:512 | Length:512 | Length:512 | Length:512 | Min. : 0.0010 | Min. : 0.0010 | Length:512 | Min. : 0.500 | Min. :-6.0350 | Min. :0.02396 | |
| Class :character | Class :character | 1st Qu.:2014 | Class :character | Class :character | Class :character | Class :character | Class :character | 1st Qu.:14.45 | Class :character | 1st Qu.: 8.000 | Class :character | Class :character | Class :character | Class :character | Class :character | 1st Qu.:100.0 | 1st Qu.: 600.0 | Class :character | Class :character | Class :character | 1st Qu.: 11 | Class :character | Class :character | 1st Qu.: 5.00 | 1st Qu.:1.000 | Class :character | 1st Qu.: 0.160 | Class :character | 1st Qu.: 0.0010 | Class :character | Class :character | 1st Qu.: 5.00 | 1st Qu.:1.000 | 1st Qu.: 0.0940 | Class :character | 1st Qu.: 0.001 | Class :character | 1st Qu.:1.000 | Class :character | Class :character | Class :character | Class :character | Class :character | Class :character | Class :character | Class :character | 1st Qu.: 0.0354 | 1st Qu.: 0.0585 | Class :character | 1st Qu.: 2.500 | 1st Qu.:-0.8778 | 1st Qu.:0.14375 | |
| Mode :character | Mode :character | Median :2019 | Mode :character | Mode :character | Mode :character | Mode :character | Mode :character | Median :18.35 | Mode :character | Median : 8.000 | Mode :character | Mode :character | Mode :character | Mode :character | Mode :character | Median :160.0 | Median : 600.0 | Mode :character | Mode :character | Mode :character | Median : 300 | Mode :character | Mode :character | Median :10.00 | Median :1.000 | Mode :character | Median : 0.298 | Mode :character | Median : 0.0100 | Mode :character | Mode :character | Median :10.00 | Median :1.000 | Median : 0.2285 | Mode :character | Median : 0.020 | Mode :character | Median :3.000 | Mode :character | Mode :character | Mode :character | Mode :character | Mode :character | Mode :character | Mode :character | Mode :character | Median : 0.1580 | Median : 0.1461 | Mode :character | Median : 5.000 | Median :-0.1671 | Median :0.14375 | |
| NA | NA | Mean :2017 | NA | NA | NA | NA | NA | Mean :21.04 | NA | Mean : 8.994 | NA | NA | NA | NA | NA | Mean :161.3 | Mean : 733.3 | NA | NA | NA | Mean : 2304 | NA | NA | Mean :11.49 | Mean :2.371 | NA | Mean : 3.494 | NA | Mean : 1.7676 | NA | NA | Mean :11.49 | Mean :2.436 | Mean : 3.2321 | NA | Mean : 1.822 | NA | Mean :2.481 | NA | NA | NA | NA | NA | NA | NA | NA | Mean : 4.4069 | Mean : 4.4491 | NA | Mean : 5.744 | Mean :-0.3631 | Mean :0.20045 | |
| NA | NA | 3rd Qu.:2019 | NA | NA | NA | NA | NA | 3rd Qu.:21.31 | NA | 3rd Qu.:11.000 | NA | NA | NA | NA | NA | 3rd Qu.:175.0 | 3rd Qu.: 900.0 | NA | NA | NA | 3rd Qu.: 300 | NA | NA | 3rd Qu.:10.00 | 3rd Qu.:5.000 | NA | 3rd Qu.: 1.083 | NA | 3rd Qu.: 0.1185 | NA | NA | 3rd Qu.:10.00 | 3rd Qu.:5.000 | 3rd Qu.: 1.0505 | NA | 3rd Qu.: 0.130 | NA | 3rd Qu.:3.000 | NA | NA | NA | NA | NA | NA | NA | NA | 3rd Qu.: 0.5600 | 3rd Qu.: 0.6516 | NA | 3rd Qu.: 5.000 | 3rd Qu.: 0.1849 | 3rd Qu.:0.28750 | |
| NA | NA | Max. :2020 | NA | NA | NA | NA | NA | Max. :79.11 | NA | Max. :14.000 | NA | NA | NA | NA | NA | Max. :300.0 | Max. :1500.0 | NA | NA | NA | Max. :59777 | NA | NA | Max. :60.00 | Max. :6.000 | NA | Max. :86.689 | NA | Max. :133.7000 | NA | NA | Max. :60.00 | Max. :6.000 | Max. :134.4379 | NA | Max. :130.500 | NA | Max. :4.000 | NA | NA | NA | NA | NA | NA | NA | NA | Max. :133.7000 | Max. :130.5000 | NA | Max. :30.000 | Max. : 3.4622 | Max. :1.43748 | |
| NA | NA | NA | NA | NA | NA | NA | NA | NA’s :284 | NA | NA | NA | NA | NA | NA | NA | NA’s :6 | NA’s :56 | NA | NA | NA | NA’s :73 | NA | NA | NA | NA | NA | NA | NA | NA’s :53 | NA | NA | NA | NA | NA | NA | NA’s :55 | NA | NA’s :198 | NA | NA | NA | NA | NA | NA | NA | NA | NA’s :330 | NA’s :328 | NA | NA | NA | NA |
Cohort_ID explains virtually no variance in the model. Hence, it was removed from the model. All the other random effects explained significant variance and were kept in subsequent models
MA_all_rand_effects <- rma.mv(lnRR, VCV_lnRR, # Add `VCV_lnRR` to account for correlated errors errors between cohorts (shared_controls)
random = list(~1|Study_ID, # Identity of the study
~1|Phylogeny, # Phylogenetic correlation
~1|Cohort_ID, # Identity of the cohort (shared controls)
~1|Species_common, # Non-phylogenetic correlation between species
~1|PFAS_type, # Type of PFAS
~1|Effect_ID), # Effect size identity
R= list(Phylogeny = cor_tree), # Assign the 'Phylogeny' argument to the phylogenetic variance-covariance matrix
test = "t",
data = dat)
summary(MA_all_rand_effects) # Cohort ID does not explain any variance ##
## Multivariate Meta-Analysis Model (k = 512; method: REML)
##
## logLik Deviance AIC BIC AICc
## -634.9177 1269.8353 1283.8353 1313.4899 1284.0580
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.5844 0.7645 10 no Study_ID no
## sigma^2.2 0.0092 0.0960 38 no Phylogeny yes
## sigma^2.3 0.0000 0.0004 153 no Cohort_ID no
## sigma^2.4 0.2081 0.4562 39 no Species_common no
## sigma^2.5 0.1009 0.3177 18 no PFAS_type no
## sigma^2.6 0.4877 0.6984 512 no Effect_ID no
##
## Test for Heterogeneity:
## Q(df = 511) = 7278.2801, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## -0.3142 0.2917 -1.0770 0.2820 -0.8874 0.2589
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MA_model <- rma.mv(lnRR, VCV_lnRR,
random = list(~1|Study_ID,
~1|Phylogeny, # Removed Cohort_ID
~1|Species_common,
~1|PFAS_type,
~1|Effect_ID),
R= list(Phylogeny = cor_tree),
test = "t",
data = dat)
summary(MA_model)##
## Multivariate Meta-Analysis Model (k = 512; method: REML)
##
## logLik Deviance AIC BIC AICc
## -634.9176 1269.8353 1281.8353 1307.2535 1282.0020
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.5844 0.7645 10 no Study_ID no
## sigma^2.2 0.0092 0.0959 38 no Phylogeny yes
## sigma^2.3 0.2081 0.4562 39 no Species_common no
## sigma^2.4 0.1009 0.3177 18 no PFAS_type no
## sigma^2.5 0.4877 0.6984 512 no Effect_ID no
##
## Test for Heterogeneity:
## Q(df = 511) = 7278.2801, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## -0.3142 0.2917 -1.0771 0.2819 -0.8874 0.2589
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## I2_total I2_Study_ID I2_Phylogeny I2_Species_common
## 0.917318637 0.385600355 0.006068127 0.137303674
## I2_PFAS_type I2_Effect_ID
## 0.066572256 0.321774224
# plot
orchard_plot(MA_model, mod = "Int", xlab = "lnRR", alpha=0.4) + # Orchard plot
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5)+ # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0, show.legend = FALSE, size = 2)+ # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_colour_manual(values = "darkorange")+ # change colours
scale_fill_manual(values="darkorange")+
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) run_model<-function(data,formula){
data<-as.data.frame(data) # convert data set into a data frame to calculate VCV matrix
VCV<-make_VCV_matrix(data, V = "var_lnRR", cluster = "Cohort_ID", obs = "Effect_ID", rho = 0.5) # create VCV matrix for the specified data
rma.mv(lnRR, VCV, # run the model, as described earlier
mods=formula,
random = list(~1|Study_ID,
~1|Phylogeny,
~1|Species_common,
~1|PFAS_type,
~1|Effect_ID),
R= list(Phylogeny = cor_tree),
test = "t",
data = data)
}plot_continuous<-function(data, model, moderator, xlab){
pred<-predict.rma(model)
data %>% mutate(fit=pred$pred,
ci.lb=pred$ci.lb,
ci.ub=pred$ci.ub,
pr.lb=pred$cr.lb,
pr.ub=pred$cr.ub) %>%
ggplot(aes(x = moderator, y = lnRR)) +
geom_ribbon(aes(ymin = pr.lb, ymax = pr.ub, color = NULL), alpha = .075) +
geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = .2) +
geom_point(aes(size=(1/sqrt(var_lnRR)), fill=Cooking_Category), shape=21, alpha=0.8) +
scale_fill_manual(values=c("#55C667FF", "goldenrod2", "dodgerblue3"))+
geom_line(aes(y = fit), size = 1.5)+
labs(x = xlab, y = "lnRR", size = "Precison (1/SE)") +
theme_bw() +
scale_size_continuous(range=c(1,9))+
geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+ # horizontal line at lnRR = 0
theme(text = element_text(size = 18, colour = "black", hjust = 0.5), # change font sizes and legend position
legend.text=element_text(size=14),
legend.position=c(0,0),
legend.justification = c(0,0),
legend.background = element_blank(),
legend.direction="horizontal",
legend.title = element_text(size=15),
panel.border=element_rect(colour="black", fill=NA, size=1.2))
}All continuous variables were z-transformed
##
## Multivariate Meta-Analysis Model (k = 512; method: REML)
##
## logLik Deviance AIC BIC AICc
## -631.9971 1263.9942 1279.9942 1313.8537 1280.2822
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.5905 0.7685 10 no Study_ID no
## sigma^2.2 0.0023 0.0481 38 no Phylogeny yes
## sigma^2.3 0.2116 0.4600 39 no Species_common no
## sigma^2.4 0.1023 0.3198 18 no PFAS_type no
## sigma^2.5 0.4881 0.6986 512 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 509) = 7233.4707, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 509) = 1.0533, p-val = 0.3686
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## Cooking_CategoryNo liquid -0.2018 0.3125 -0.6457 0.5187 -0.8159 0.4122
## Cooking_Categoryoil-based -0.3712 0.2971 -1.2495 0.2121 -0.9548 0.2124
## Cooking_Categorywater-based -0.2932 0.2950 -0.9939 0.3207 -0.8729 0.2864
##
## Cooking_CategoryNo liquid
## Cooking_Categoryoil-based
## Cooking_Categorywater-based
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R2_marginal R2_coditional
## 0.002563516 0.650990549
# plot
orchard_plot(category_model, mod = "Cooking_Category", xlab = "lnRR", alpha=0.4)+
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5)+ # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0, show.legend = FALSE, size = 2)+ # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_colour_manual(values = c("#55C667FF", "goldenrod2", "dodgerblue3"))+ # change colours
scale_fill_manual(values=c("#55C667FF", "goldenrod2", "dodgerblue3"))+
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13))##
## Multivariate Meta-Analysis Model (k = 512; method: REML)
##
## logLik Deviance AIC BIC AICc
## -633.4258 1266.8515 1280.8515 1310.4924 1281.0746
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.5830 0.7636 10 no Study_ID no
## sigma^2.2 0.0106 0.1028 38 no Phylogeny yes
## sigma^2.3 0.2085 0.4566 39 no Species_common no
## sigma^2.4 0.1061 0.3257 18 no PFAS_type no
## sigma^2.5 0.4880 0.6985 512 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 510) = 7223.9798, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 510) = 0.1431, p-val = 0.7054
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -0.4218 0.4087 -1.0322 0.3024 -1.2247 0.3810
## PFAS_carbon_chain 0.0119 0.0315 0.3783 0.7054 -0.0499 0.0738
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R2_marginal R2_coditional
## 0.0005515213 0.6506803497
# Temperature_in_Celsius
temp_model<-run_model(dat, ~scale(Temperature_in_Celsius)) # z-transformed
summary(temp_model)##
## Multivariate Meta-Analysis Model (k = 506; method: REML)
##
## logLik Deviance AIC BIC AICc
## -626.2464 1252.4927 1266.4927 1296.0508 1266.7185
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.5805 0.7619 10 no Study_ID no
## sigma^2.2 0.0054 0.0735 38 no Phylogeny yes
## sigma^2.3 0.2079 0.4559 39 no Species_common no
## sigma^2.4 0.0974 0.3122 18 no PFAS_type no
## sigma^2.5 0.4896 0.6997 506 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 504) = 7121.6638, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 504) = 0.0318, p-val = 0.8586
##
## Model Results:
##
## estimate se tval pval ci.lb
## intrcpt -0.3001 0.2911 -1.0309 0.3031 -0.8719
## scale(Temperature_in_Celsius) 0.0144 0.0810 0.1783 0.8586 -0.1447
## ci.ub
## intrcpt 0.2718
## scale(Temperature_in_Celsius) 0.1736
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R2_marginal R2_coditional
## 0.000150956 0.645470163
# Plot
dat.temp<-filter(dat, Temperature_in_Celsius!="NA")
plot_continuous(dat.temp, temp_model, dat.temp$Temperature_in_Celsius, "Cooking temperature")# Length_cooking_time_in_s
time_model<-run_model(dat, ~scale(Length_cooking_time_in_s)) # z-transformed
summary(time_model)##
## Multivariate Meta-Analysis Model (k = 456; method: REML)
##
## logLik Deviance AIC BIC AICc
## -525.7357 1051.4714 1065.4714 1094.2980 1065.7225
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.5387 0.7340 9 no Study_ID no
## sigma^2.2 0.0000 0.0001 30 no Phylogeny yes
## sigma^2.3 0.1425 0.3775 30 no Species_common no
## sigma^2.4 0.1009 0.3176 17 no PFAS_type no
## sigma^2.5 0.3964 0.6296 456 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 454) = 3999.2874, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 454) = 24.7942, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb
## intrcpt -0.5531 0.2884 -1.9175 0.0558 -1.1199
## scale(Length_cooking_time_in_s) -0.2646 0.0531 -4.9794 <.0001 -0.3690
## ci.ub
## intrcpt 0.0138 .
## scale(Length_cooking_time_in_s) -0.1601 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R2_marginal R2_coditional
## 0.05605974 0.68251140
# Plot
dat.time<-filter(dat, Length_cooking_time_in_s!="NA")
plot_continuous(dat.time, time_model, dat.time$Length_cooking_time_in_s, "Cooking time (s)")# Volume_liquid_ml
volume_model<-run_model(dat, ~scale(log(Volume_liquid_ml))) # logged and z-transformed
summary(volume_model)##
## Multivariate Meta-Analysis Model (k = 439; method: REML)
##
## logLik Deviance AIC BIC AICc
## -552.0542 1104.1084 1118.1084 1146.6680 1118.3695
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.5079 0.7126 8 no Study_ID no
## sigma^2.2 0.0048 0.0692 34 no Phylogeny yes
## sigma^2.3 0.1498 0.3870 35 no Species_common no
## sigma^2.4 0.1177 0.3431 18 no PFAS_type no
## sigma^2.5 0.5153 0.7178 439 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 437) = 5805.2399, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 437) = 5.9117, p-val = 0.0154
##
## Model Results:
##
## estimate se tval pval ci.lb
## intrcpt -0.3568 0.2978 -1.1981 0.2315 -0.9421
## scale(log(Volume_liquid_ml)) -0.2543 0.1046 -2.4314 0.0154 -0.4599
## ci.ub
## intrcpt 0.2285
## scale(log(Volume_liquid_ml)) -0.0487 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R2_marginal R2_coditional
## 0.0475590 0.6211531
# Plot
dat.volume<-filter(dat, Volume_liquid_ml!="NA")
plot_continuous(dat.volume, volume_model, log(dat.volume$Volume_liquid_ml), "Volume of liquid (mL)")# Moisture_loss_in_percent
moisture_model<-run_model(dat, ~scale(Moisture_loss_in_percent))
summary(moisture_model)##
## Multivariate Meta-Analysis Model (k = 228; method: REML)
##
## logLik Deviance AIC BIC AICc
## -234.1714 468.3428 482.3428 506.2865 482.8566
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0775 0.2785 6 no Study_ID no
## sigma^2.2 0.2316 0.4812 18 no Phylogeny yes
## sigma^2.3 0.1307 0.3615 18 no Species_common no
## sigma^2.4 0.0094 0.0968 17 no PFAS_type no
## sigma^2.5 0.3220 0.5674 228 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 226) = 2492.6080, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 226) = 0.0295, p-val = 0.8638
##
## Model Results:
##
## estimate se tval pval ci.lb
## intrcpt 0.5347 0.3311 1.6147 0.1078 -0.1178
## scale(Moisture_loss_in_percent) -0.0117 0.0683 -0.1717 0.8638 -0.1463
## ci.ub
## intrcpt 1.1872
## scale(Moisture_loss_in_percent) 0.1229
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R2_marginal R2_coditional
## 0.0001783538 0.5825498843
# Plot
dat.moisture<-filter(dat, Moisture_loss_in_percent!="NA")
plot_continuous(dat.moisture, moisture_model, dat.moisture$Moisture_loss_in_percent, "Percentage of moisture loss")# Full_model
full_model <- run_model(dat, ~ - 1 +
Cooking_Category +
scale(Temperature_in_Celsius) +
scale(Length_cooking_time_in_s) +
scale(PFAS_carbon_chain) +
scale(log(Volume_liquid_ml)))
summary(full_model)##
## Multivariate Meta-Analysis Model (k = 399; method: REML)
##
## logLik Deviance AIC BIC AICc
## -442.1936 884.3873 908.3873 956.0424 909.2105
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0835 0.2890 7 no Study_ID no
## sigma^2.2 0.0000 0.0001 26 no Phylogeny yes
## sigma^2.3 0.0941 0.3067 26 no Species_common no
## sigma^2.4 0.1212 0.3481 17 no PFAS_type no
## sigma^2.5 0.3829 0.6188 399 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 392) = 2935.7864, p-val < .0001
##
## Test of Moderators (coefficients 1:7):
## F(df1 = 7, df2 = 392) = 13.2854, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb
## Cooking_CategoryNo liquid -0.5772 0.2836 -2.0356 0.0425 -1.1348
## Cooking_Categoryoil-based -0.7378 0.1891 -3.9010 0.0001 -1.1097
## Cooking_Categorywater-based -0.4054 0.2202 -1.8409 0.0664 -0.8384
## scale(Temperature_in_Celsius) -0.0147 0.1119 -0.1317 0.8953 -0.2348
## scale(Length_cooking_time_in_s) -0.4005 0.0577 -6.9370 <.0001 -0.5140
## scale(PFAS_carbon_chain) 0.0619 0.0799 0.7746 0.4390 -0.0952
## scale(log(Volume_liquid_ml)) -0.7214 0.1027 -7.0271 <.0001 -0.9233
## ci.ub
## Cooking_CategoryNo liquid -0.0197 *
## Cooking_Categoryoil-based -0.3660 ***
## Cooking_Categorywater-based 0.0276 .
## scale(Temperature_in_Celsius) 0.2053
## scale(Length_cooking_time_in_s) -0.2870 ***
## scale(PFAS_carbon_chain) 0.2189
## scale(log(Volume_liquid_ml)) -0.5196 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R2_marginal R2_coditional
## 0.4600194 0.6967281
## Cooking_CategoryNo liquid Cooking_Categoryoil-based
## 1.6588 3.3203
## Cooking_Categorywater-based scale(Temperature_in_Celsius)
## 4.1165 2.3290
## scale(Length_cooking_time_in_s) scale(PFAS_carbon_chain)
## 1.0993 1.0018
## scale(log(Volume_liquid_ml))
## 1.3527
dat %>% select(Temperature_in_Celsius, Length_cooking_time_in_s, PFAS_carbon_chain, Volume_liquid_ml) %>%
ggpairs()SHOULD WE USE PROPORTIONAL OR EQUAL WEIGHTS????
Inspection of the plots highlighted potential significant decreases in PFAS content with increased cooking time and volume of cooking. Hence, here we used emmeans (download from remotes::install_github(“rvlenth/emmeans”, dependencies = TRUE, build_opts = "")) to generate marginalised means at specified values of the different predictors. Such analysis enable the quantification of the mean effect size after controlling for different values of the moderators.
# Full model in original units (not z-transformation)
dat$log_Volume_liquid_ml<-log(dat$Volume_liquid_ml)
full_model_org_units <- run_model(dat, ~ - 1 +
Cooking_Category +
Temperature_in_Celsius +
Length_cooking_time_in_s +
PFAS_carbon_chain +
log_Volume_liquid_ml)
save(full_model, file = here("Rdata", "full_model_org_units.RData"))res_mean_values <- qdrg(object = full_model_org_units, data = dat) # Generate an grid of marginal means at the mean of each predictor variable
res_mean_values## 'emmGrid' object with variables:
## Cooking_Category = No liquid, oil-based, water-based
## Temperature_in_Celsius = 162.42
## Length_cooking_time_in_s = 726.77
## PFAS_carbon_chain = 9.0226
## log_Volume_liquid_ml = 4.8171
## 1 emmean SE df asymp.LCL asymp.UCL
## overall -0.556 0.196 Inf -0.939 -0.172
##
## Results are averaged over the levels of: Cooking_Category
## Confidence level used: 0.95
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid -0.559 0.283 Inf -1.115 -0.00371
## oil-based -0.720 0.189 Inf -1.090 -0.34979
## water-based -0.387 0.221 Inf -0.821 0.04647
##
## Confidence level used: 0.95
Here, we generate estimates at cooking times of 2, 5, 10, 15, 20 and 25 min.
res_cooking_time <- qdrg(object = full_model_org_units, data = dat, at = list(Length_cooking_time_in_s = c(120, 300, 600, 900, 1200, 1500)))
res_cooking_time ## 'emmGrid' object with variables:
## Cooking_Category = No liquid, oil-based, water-based
## Temperature_in_Celsius = 162.42
## Length_cooking_time_in_s = 120, 300, 600, 900, 1200, 1500
## PFAS_carbon_chain = 9.0226
## log_Volume_liquid_ml = 4.8171
## Length_cooking_time_in_s = 120:
## emmean SE df asymp.LCL asymp.UCL
## 0.2760 0.225 Inf -0.165 0.71677
##
## Length_cooking_time_in_s = 300:
## emmean SE df asymp.LCL asymp.UCL
## 0.0293 0.210 Inf -0.381 0.44013
##
## Length_cooking_time_in_s = 600:
## emmean SE df asymp.LCL asymp.UCL
## -0.3818 0.196 Inf -0.766 0.00268
##
## Length_cooking_time_in_s = 900:
## emmean SE df asymp.LCL asymp.UCL
## -0.7930 0.200 Inf -1.185 -0.40066
##
## Length_cooking_time_in_s = 1200:
## emmean SE df asymp.LCL asymp.UCL
## -1.2041 0.221 Inf -1.636 -0.77173
##
## Length_cooking_time_in_s = 1500:
## emmean SE df asymp.LCL asymp.UCL
## -1.6152 0.254 Inf -2.112 -1.11827
##
## Results are averaged over the levels of: Cooking_Category
## Confidence level used: 0.95
## Length_cooking_time_in_s = 120:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid 0.2723 0.305 Inf -0.3262 0.871
## oil-based 0.1117 0.217 Inf -0.3140 0.537
## water-based 0.4441 0.248 Inf -0.0419 0.930
##
## Length_cooking_time_in_s = 300:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid 0.0256 0.294 Inf -0.5506 0.602
## oil-based -0.1350 0.202 Inf -0.5307 0.261
## water-based 0.1974 0.234 Inf -0.2614 0.656
##
## Length_cooking_time_in_s = 600:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid -0.3856 0.284 Inf -0.9422 0.171
## oil-based -0.5462 0.189 Inf -0.9164 -0.176
## water-based -0.2137 0.222 Inf -0.6487 0.221
##
## Length_cooking_time_in_s = 900:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid -0.7967 0.286 Inf -1.3578 -0.236
## oil-based -0.9573 0.194 Inf -1.3376 -0.577
## water-based -0.6249 0.225 Inf -1.0664 -0.183
##
## Length_cooking_time_in_s = 1200:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid -1.2078 0.300 Inf -1.7967 -0.619
## oil-based -1.3684 0.216 Inf -1.7917 -0.945
## water-based -1.0360 0.243 Inf -1.5131 -0.559
##
## Length_cooking_time_in_s = 1500:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid -1.6190 0.325 Inf -2.2558 -0.982
## oil-based -1.7796 0.250 Inf -2.2701 -1.289
## water-based -1.4472 0.273 Inf -1.9832 -0.911
##
## Confidence level used: 0.95
Here, we generate marginalised estimates at volumes of liquid of ~10, 500, 1000, 5000 and 10000 mL. We did not look at the means for different cooking categories because they are inherently different in the volume of liquid used.
res_volume <- qdrg(object = full_model_org_units, data = dat, at = list(log_Volume_liquid_ml= c(2.3, 6.2, 6.9, 8.5, 9.2)))
res_volume ## 'emmGrid' object with variables:
## Cooking_Category = No liquid, oil-based, water-based
## Temperature_in_Celsius = 162.42
## Length_cooking_time_in_s = 726.77
## PFAS_carbon_chain = 9.0226
## log_Volume_liquid_ml = 2.3, 6.2, 6.9, 8.5, 9.2
## log_Volume_liquid_ml = 2.3:
## emmean SE df asymp.LCL asymp.UCL
## 0.313 0.247 Inf -0.171 0.797
##
## log_Volume_liquid_ml = 6.2:
## emmean SE df asymp.LCL asymp.UCL
## -1.033 0.197 Inf -1.419 -0.646
##
## log_Volume_liquid_ml = 6.9:
## emmean SE df asymp.LCL asymp.UCL
## -1.274 0.207 Inf -1.679 -0.869
##
## log_Volume_liquid_ml = 8.5:
## emmean SE df asymp.LCL asymp.UCL
## -1.826 0.245 Inf -2.307 -1.345
##
## log_Volume_liquid_ml = 9.2:
## emmean SE df asymp.LCL asymp.UCL
## -2.068 0.268 Inf -2.593 -1.542
##
## Results are averaged over the levels of: Cooking_Category
## Confidence level used: 0.95
Here, we generate marginalised estimates for PFAS of 3, 6, 9 and 12 carbon chains
res_PFAS <- qdrg(object = full_model_org_units, data = dat, at = list(PFAS_carbon_chain= c(3, 6, 9, 12)))
res_PFAS ## 'emmGrid' object with variables:
## Cooking_Category = No liquid, oil-based, water-based
## Temperature_in_Celsius = 162.42
## Length_cooking_time_in_s = 726.77
## PFAS_carbon_chain = 3, 6, 9, 12
## log_Volume_liquid_ml = 4.8171
## PFAS_carbon_chain = 3:
## emmean SE df asymp.LCL asymp.UCL
## -0.715 0.288 Inf -1.281 -0.1502
##
## PFAS_carbon_chain = 6:
## emmean SE df asymp.LCL asymp.UCL
## -0.636 0.224 Inf -1.075 -0.1968
##
## PFAS_carbon_chain = 9:
## emmean SE df asymp.LCL asymp.UCL
## -0.556 0.196 Inf -0.940 -0.1726
##
## PFAS_carbon_chain = 12:
## emmean SE df asymp.LCL asymp.UCL
## -0.476 0.218 Inf -0.904 -0.0489
##
## Results are averaged over the levels of: Cooking_Category
## Confidence level used: 0.95
## PFAS_carbon_chain = 3:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid -0.719 0.354 Inf -1.412 -0.02615
## oil-based -0.880 0.283 Inf -1.435 -0.32437
## water-based -0.547 0.307 Inf -1.149 0.05461
##
## PFAS_carbon_chain = 6:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid -0.640 0.304 Inf -1.234 -0.04466
## oil-based -0.800 0.218 Inf -1.227 -0.37337
## water-based -0.468 0.247 Inf -0.952 0.01678
##
## PFAS_carbon_chain = 9:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid -0.560 0.283 Inf -1.116 -0.00428
## oil-based -0.720 0.189 Inf -1.091 -0.35035
## water-based -0.388 0.221 Inf -0.822 0.04591
##
## PFAS_carbon_chain = 12:
## Cooking_Category emmean SE df asymp.LCL asymp.UCL
## No liquid -0.480 0.300 Inf -1.067 0.10694
## oil-based -0.641 0.212 Inf -1.057 -0.22475
## water-based -0.308 0.241 Inf -0.781 0.16401
##
## Confidence level used: 0.95
Here, we investigated whether the effect of the continuous moderators on lnRR vary depending on the cooking category. Hence, we performed subset analyses for each cooking category.
oil_dat<-filter(dat, Cooking_Category=="oil-based")
include <- row.names(cor_tree) %in% oil_dat$Phylogeny # Check which rows are present in the phylogenetic tree
cor_tree_oil <- cor_tree[include, include] # Only include the species that match the reduced data set
run_model_oil<-function(data,formula){
data<-as.data.frame(data) # convert data set into a data frame to calculate VCV matrix
VCV<-make_VCV_matrix(data, V = "var_lnRR", cluster = "Cohort_ID", obs = "Effect_ID", rho = 0.5) # create VCV matrix for the specified data
rma.mv(lnRR, VCV, # run the model, as described earlier
mods=formula,
random = list(~1|Study_ID,
~1|Phylogeny,
~1|Species_common,
~1|PFAS_type,
~1|Effect_ID),
R= list(Phylogeny = cor_tree_oil), # cor_tree_oil here
test = "t",
data = data)
}full_model_oil <- run_model_oil(oil_dat, ~ scale(Temperature_in_Celsius) +
scale(Length_cooking_time_in_s) +
scale(PFAS_carbon_chain) +
scale(log(Volume_liquid_ml)))
summary(full_model_oil)##
## Multivariate Meta-Analysis Model (k = 263; method: REML)
##
## logLik Deviance AIC BIC AICc
## -184.0059 368.0118 388.0118 423.5414 388.9025
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0000 0.0001 6 no Study_ID no
## sigma^2.2 0.0000 0.0000 19 no Phylogeny yes
## sigma^2.3 0.0141 0.1188 19 no Species_common no
## sigma^2.4 0.0433 0.2080 16 no PFAS_type no
## sigma^2.5 0.1124 0.3353 263 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 258) = 573.2766, p-val < .0001
##
## Test of Moderators (coefficients 2:5):
## F(df1 = 4, df2 = 258) = 27.1829, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb
## intrcpt -0.4370 0.0959 -4.5554 <.0001 -0.6259
## scale(Temperature_in_Celsius) -0.0039 0.0817 -0.0474 0.9622 -0.1647
## scale(Length_cooking_time_in_s) -0.3805 0.0485 -7.8388 <.0001 -0.4761
## scale(PFAS_carbon_chain) 0.1286 0.0613 2.0957 0.0371 0.0078
## scale(log(Volume_liquid_ml)) -0.4032 0.0893 -4.5131 <.0001 -0.5791
## ci.ub
## intrcpt -0.2481 ***
## scale(Temperature_in_Celsius) 0.1570
## scale(Length_cooking_time_in_s) -0.2849 ***
## scale(PFAS_carbon_chain) 0.2494 *
## scale(log(Volume_liquid_ml)) -0.2273 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
water_dat<-filter(dat, Cooking_Category=="water-based")
include <- row.names(cor_tree) %in% water_dat$Phylogeny # Check which rows are present in the phylogenetic tree
cor_tree_water <- cor_tree[include, include] # Only include the species that match the reduced data set
run_model_water<-function(data,formula){
data<-as.data.frame(data) # convert data set into a data frame to calculate VCV matrix
VCV<-make_VCV_matrix(data, V = "var_lnRR", cluster = "Cohort_ID", obs = "Effect_ID", rho = 0.5) # create VCV matrix for the specified data
rma.mv(lnRR, VCV, # run the model, as described earlier
mods=formula,
random = list(~1|Study_ID,
~1|Phylogeny,
~1|Species_common,
~1|PFAS_type,
~1|Effect_ID),
R= list(Phylogeny = cor_tree_water), # cor_tree_water here
test = "t",
data = data)
}full_model_water <- run_model_water(water_dat, ~ scale(Temperature_in_Celsius) +
scale(Length_cooking_time_in_s) +
scale(PFAS_carbon_chain) +
scale(log(Volume_liquid_ml)))
summary(full_model_water)##
## Multivariate Meta-Analysis Model (k = 121; method: REML)
##
## logLik Deviance AIC BIC AICc
## -179.5070 359.0139 377.0139 401.8735 378.6961
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.3578 0.5982 6 no Study_ID no
## sigma^2.2 0.0000 0.0002 19 no Phylogeny yes
## sigma^2.3 0.0000 0.0039 19 no Species_common no
## sigma^2.4 0.4043 0.6359 15 no PFAS_type no
## sigma^2.5 0.9470 0.9732 121 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 117) = 2237.4353, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 117) = 4.6171, p-val = 0.0043
##
## Model Results:
##
## estimate se tval pval ci.lb
## intrcpt -0.9408 0.3396 -2.7700 0.0065 -1.6134
## scale(Length_cooking_time_in_s) -0.4503 0.1591 -2.8307 0.0055 -0.7653
## scale(PFAS_carbon_chain) -0.0082 0.1688 -0.0487 0.9612 -0.3426
## scale(log(Volume_liquid_ml)) -0.7986 0.2779 -2.8732 0.0048 -1.3491
## ci.ub
## intrcpt -0.2682 **
## scale(Length_cooking_time_in_s) -0.1353 **
## scale(PFAS_carbon_chain) 0.3261
## scale(log(Volume_liquid_ml)) -0.2481 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Not very relevant because all effect sizes are from one study here. Also, the model does not converge when adding VCV_lnRR
dry_dat<-filter(dat, Cooking_Category=="No liquid")
include <- row.names(cor_tree) %in% dry_dat$Phylogeny # Check which rows are present in the phylogenetic tree
cor_tree_dry <- cor_tree[include, include] # Only include the species that match the reduced data set
run_model_dry<-function(data,formula){
data<-as.data.frame(data) # convert data set into a data frame to calculate VCV matrix
rma.mv(lnRR, var_lnRR, # run the model with var_lnRR instead of lnCVR
mods=formula,
random = list(~1|Study_ID,
~1|Phylogeny,
~1|Species_common,
~1|PFAS_type,
~1|Effect_ID),
R= list(Phylogeny = cor_tree_dry), # cor_tree_dry here
test = "t",
data = data)
}full_model_dry <- run_model_dry(dry_dat, ~ scale(Length_cooking_time_in_s)) # Nodel does not converge with VCV_lnRR
summary(full_model_dry)##
## Multivariate Meta-Analysis Model (k = 47; method: REML)
##
## logLik Deviance AIC BIC AICc
## -12.9722 25.9445 37.9445 48.7844 40.1550
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0000 0.0000 1 yes Study_ID no
## sigma^2.2 0.0046 0.0679 8 no Phylogeny yes
## sigma^2.3 0.0022 0.0471 8 no Species_common no
## sigma^2.4 0.0735 0.2711 2 no PFAS_type no
## sigma^2.5 0.0000 0.0000 47 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 45) = 40.1184, p-val = 0.6785
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 45) = 38.2787, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb
## intrcpt -0.7770 0.2071 -3.7513 0.0005 -1.1942
## scale(Length_cooking_time_in_s) -0.3461 0.0559 -6.1870 <.0001 -0.4588
## ci.ub
## intrcpt -0.3598 ***
## scale(Length_cooking_time_in_s) -0.2334 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
oil_dat <- filter(dat, Cooking_Category=="oil-based")
water_dat <- filter(dat, Cooking_Category=="water-based")
dry_dat <- filter(dat, Cooking_Category=="No liquid")
oil_dat_time<-filter(oil_dat, Length_cooking_time_in_s!="NA")
water_dat_time<-filter(water_dat, Length_cooking_time_in_s!="NA")
dry_dat_time<-filter(dry_dat, Length_cooking_time_in_s!="NA")
model_oil_time<-run_model_oil(oil_dat_time, ~Length_cooking_time_in_s)
model_water_time<-run_model_water(water_dat_time, ~Length_cooking_time_in_s)
model_dry_time<-run_model_dry(dry_dat_time, ~Length_cooking_time_in_s)
pred_oil_time<-predict.rma(model_oil_time)
pred_water_time<-predict.rma(model_water_time)
pred_dry_time<-predict.rma(model_dry_time)
oil_dat_time<-mutate(oil_dat_time,
ci.lb = pred_oil_time$ci.lb, # lower bound of the confidence interval for oil
ci.ub = pred_oil_time$ci.ub, # upper bound of the confidence interval for oil
fit = pred_oil_time$pred) # regression line for oil
water_dat_time<-mutate(water_dat_time,
ci.lb = pred_water_time$ci.lb, # lower bound of the confidence interval for water
ci.ub = pred_water_time$ci.ub, # upper bound of the confidence interval for water
fit = pred_water_time$pred) # regression line for water
dry_dat_time<-mutate(dry_dat_time,
ci.lb = pred_dry_time$ci.lb, # lower bound of the confidence interval for dry
ci.ub = pred_dry_time$ci.ub, # upper bound of the confidence interval for dry
fit = pred_dry_time$pred) # regression line for dryggplot(dat,aes(x = Length_cooking_time_in_s, y = lnRR, fill=Cooking_Category)) +
geom_ribbon(data=water_dat_time, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=water_dat_time,aes(y = fit), size = 1.5, col="dodgerblue")+
geom_ribbon(data=oil_dat_time, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.3) +
geom_line(data=oil_dat_time,aes(y = fit), size = 1.5, col="goldenrod")+
geom_ribbon(data=dry_dat_time, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=dry_dat_time,aes(y = fit), size = 1.5, col="palegreen3")+
geom_point(aes(size=(1/sqrt(var_lnRR)), fill=Cooking_Category), shape=21, alpha=0.8) +
scale_fill_manual(values=c("#55C667FF", "goldenrod2", "dodgerblue3"))+
labs(x = "Cooking time (s)", y = "lnRR", size = "Precison (1/SE)") +
scale_size_continuous(range=c(1,10))+
theme_bw() +
geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+ # horizontal line at lnRR = 0
theme(text = element_text(size = 18, colour = "black", hjust = 0.5), # change font sizes and legend position
legend.text=element_text(size=14),
legend.position=c(0,0),
legend.justification = c(0,0),
legend.background = element_blank(),
legend.direction="horizontal",
legend.title = element_text(size=15),
panel.border=element_rect(colour="black", fill=NA, size=1.2)) oil_dat_vol<-filter(oil_dat, Volume_liquid_ml!="NA")
water_dat_vol<-filter(water_dat, Volume_liquid_ml!="NA")
dry_dat_vol<-filter(dry_dat, Volume_liquid_ml!="NA")
model_oil_vol<-run_model_oil(oil_dat_vol, ~log(Volume_liquid_ml))
model_water_vol<-run_model_water(water_dat_vol, ~log(Volume_liquid_ml))
model_dry_vol<-run_model_dry(dry_dat_vol, ~log(Volume_liquid_ml))
pred_oil_vol<-predict.rma(model_oil_vol)
pred_water_vol<-predict.rma(model_water_vol)
pred_dry_vol<-predict.rma(model_dry_vol)
oil_dat_vol<-mutate(oil_dat_vol,
ci.lb = pred_oil_vol$ci.lb,
ci.ub = pred_oil_vol$ci.ub,
fit = pred_oil_vol$pred)
water_dat_vol<-mutate(water_dat_vol,
ci.lb = pred_water_vol$ci.lb,
ci.ub = pred_water_vol$ci.ub,
fit = pred_water_vol$pred)
dry_dat_vol<-mutate(dry_dat_vol,
ci.lb = pred_dry_vol$ci.lb,
ci.ub = pred_dry_vol$ci.ub,
fit = pred_dry_vol$pred) ggplot(dat,aes(x = log(Volume_liquid_ml), y = lnRR, fill=Cooking_Category)) +
geom_ribbon(data=water_dat_vol, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=water_dat_vol,aes(y = fit), size = 1.5, col="dodgerblue")+
geom_ribbon(data=oil_dat_vol, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.3) +
geom_line(data=oil_dat_vol,aes(y = fit), size = 1.5, col="goldenrod")+
geom_ribbon(data=dry_dat_vol, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=dry_dat_vol,aes(y = fit), size = 1.5, col="palegreen3")+
geom_point(aes(size=(1/sqrt(var_lnRR)), fill=Cooking_Category), shape=21, alpha=0.8) +
scale_fill_manual(values=c("#55C667FF", "goldenrod2", "dodgerblue3"))+
labs(x = "ln(Volume of liquid (mL))", y = "lnRR", size = "Precison (1/SE)") +
scale_size_continuous(range=c(1,10))+
theme_bw() +
geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+
theme(text = element_text(size = 18, colour = "black", hjust = 0.5),
legend.text=element_text(size=14),
legend.position=c(0,0),
legend.justification = c(0,0),
legend.background = element_blank(),
legend.direction="horizontal",
legend.title = element_text(size=15),
panel.border=element_rect(colour="black", fill=NA, size=1.2)) oil_dat_PFAS<-filter(oil_dat, PFAS_carbon_chain!="NA")
water_dat_PFAS<-filter(water_dat, PFAS_carbon_chain!="NA")
dry_dat_PFAS<-filter(dry_dat, PFAS_carbon_chain!="NA")
model_oil_PFAS<-run_model_oil(oil_dat_PFAS, ~PFAS_carbon_chain)
model_water_PFAS<-run_model_water(water_dat_PFAS, ~PFAS_carbon_chain)
model_dry_PFAS<-run_model_dry(dry_dat_PFAS, ~PFAS_carbon_chain)
pred_oil_PFAS<-predict.rma(model_oil_PFAS)
pred_water_PFAS<-predict.rma(model_water_PFAS)
pred_dry_PFAS<-predict.rma(model_dry_PFAS)
oil_dat_PFAS<-mutate(oil_dat_PFAS,
ci.lb = pred_oil_PFAS$ci.lb,
ci.ub = pred_oil_PFAS$ci.ub,
fit = pred_oil_PFAS$pred)
water_dat_PFAS<-mutate(water_dat_PFAS,
ci.lb = pred_water_PFAS$ci.lb,
ci.ub = pred_water_PFAS$ci.ub,
fit = pred_water_PFAS$pred)
dry_dat_PFAS<-mutate(dry_dat_PFAS,
ci.lb = pred_dry_PFAS$ci.lb,
ci.ub = pred_dry_PFAS$ci.ub,
fit = pred_dry_PFAS$pred) ggplot(dat,aes(x = PFAS_carbon_chain, y = lnRR, fill=Cooking_Category)) +
geom_ribbon(data=dry_dat_PFAS, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=dry_dat_PFAS,aes(y = fit), size = 1.5, col="palegreen3")+
geom_ribbon(data=water_dat_PFAS, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=water_dat_PFAS,aes(y = fit), size = 1.5, col="dodgerblue")+
geom_ribbon(data=oil_dat_PFAS, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.3) +
geom_line(data=oil_dat_PFAS,aes(y = fit), size = 1.5, col="goldenrod")+
geom_point(aes(size=(1/sqrt(var_lnRR)), fill=Cooking_Category), shape=21, alpha=0.8) +
scale_fill_manual(values=c("#55C667FF", "goldenrod2", "dodgerblue3"))+
labs(x = "PFAS carbon chain length", y = "lnRR", size = "Precison (1/SE)") +
scale_size_continuous(range=c(1,10))+
theme_bw() +
geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+
theme(text = element_text(size = 18, colour = "black", hjust = 0.5),
legend.text=element_text(size=14),
legend.position=c(0,0),
legend.justification = c(0,0),
legend.background = element_blank(),
legend.direction="horizontal",
legend.title = element_text(size=15),
panel.border=element_rect(colour="black", fill=NA, size=1.2))egger_all <- run_model(dat, ~ - 1 + Cooking_Category +
I(sqrt(1/N_tilde)) +
scale(Publication_year) +
scale(Temperature_in_Celsius) +
scale(Length_cooking_time_in_s) +
scale(PFAS_carbon_chain) +
scale(log(Volume_liquid_ml)))
summary(egger_all)##
## Multivariate Meta-Analysis Model (k = 399; method: REML)
##
## logLik Deviance AIC BIC AICc
## -437.7512 875.5024 903.5024 959.0284 904.6224
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.1374 0.3707 7 no Study_ID no
## sigma^2.2 0.0000 0.0001 26 no Phylogeny yes
## sigma^2.3 0.0082 0.0907 26 no Species_common no
## sigma^2.4 0.1188 0.3447 17 no PFAS_type no
## sigma^2.5 0.3978 0.6307 399 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 390) = 2866.3281, p-val < .0001
##
## Test of Moderators (coefficients 1:9):
## F(df1 = 9, df2 = 390) = 9.9236, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb
## Cooking_CategoryNo liquid -0.2088 0.3994 -0.5227 0.6015 -0.9940
## Cooking_Categoryoil-based -0.3346 0.3457 -0.9679 0.3337 -1.0141
## Cooking_Categorywater-based 0.0285 0.3478 0.0821 0.9346 -0.6553
## I(sqrt(1/N_tilde)) -0.7252 0.5313 -1.3649 0.1731 -1.7699
## scale(Publication_year) 0.1942 0.1222 1.5890 0.1129 -0.0461
## scale(Temperature_in_Celsius) 0.0243 0.1107 0.2195 0.8263 -0.1934
## scale(Length_cooking_time_in_s) -0.3939 0.0583 -6.7509 <.0001 -0.5086
## scale(PFAS_carbon_chain) 0.0708 0.0800 0.8849 0.3767 -0.0865
## scale(log(Volume_liquid_ml)) -0.6800 0.1037 -6.5596 <.0001 -0.8839
## ci.ub
## Cooking_CategoryNo liquid 0.5765
## Cooking_Categoryoil-based 0.3450
## Cooking_Categorywater-based 0.7124
## I(sqrt(1/N_tilde)) 0.3194
## scale(Publication_year) 0.4346
## scale(Temperature_in_Celsius) 0.2420
## scale(Length_cooking_time_in_s) -0.2792 ***
## scale(PFAS_carbon_chain) 0.2282
## scale(log(Volume_liquid_ml)) -0.4762 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#funnel(egger_all, yaxis = "seinv")
# little evidence
egger_n <- run_model(dat, ~ I(sqrt(1/N_tilde)))
summary(egger_n)##
## Multivariate Meta-Analysis Model (k = 512; method: REML)
##
## logLik Deviance AIC BIC AICc
## -632.8391 1265.6782 1279.6782 1309.3191 1279.9013
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.6010 0.7752 10 no Study_ID no
## sigma^2.2 0.0044 0.0664 38 no Phylogeny yes
## sigma^2.3 0.1987 0.4458 39 no Species_common no
## sigma^2.4 0.1008 0.3175 18 no PFAS_type no
## sigma^2.5 0.4887 0.6991 512 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 510) = 6790.0424, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 510) = 0.6334, p-val = 0.4265
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -0.0930 0.4055 -0.2294 0.8186 -0.8896 0.7036
## I(sqrt(1/N_tilde)) -0.5005 0.6289 -0.7959 0.4265 -1.7361 0.7350
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multivariate Meta-Analysis Model (k = 512; method: REML)
##
## logLik Deviance AIC BIC AICc
## -631.7358 1263.4716 1277.4716 1307.1125 1277.6947
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.5567 0.7461 10 no Study_ID no
## sigma^2.2 0.0145 0.1206 38 no Phylogeny yes
## sigma^2.3 0.2046 0.4524 39 no Species_common no
## sigma^2.4 0.1014 0.3184 18 no PFAS_type no
## sigma^2.5 0.4878 0.6984 512 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 510) = 7278.1828, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 510) = 1.2853, p-val = 0.2574
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -165.8555 146.0186 -1.1359 0.2566 -452.7275 121.0165
## Publication_year 0.0821 0.0724 1.1337 0.2574 -0.0602 0.2243
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
Here, we iteratively removed one study at the time and investigated how it affects the overall mean. Removing one of the study particularly modifies the estimate, but none of these models show a significant overall difference in PFAS concentration with cooking.
dat$Study_ID<-as.factor(dat$Study_ID)
dat<-as.data.frame(dat) # Only work with a dataframe
VCV_matrix<-list() # will need new VCV matrices because the sample size will be iteratively reduced
Leave1studyout<-list() # create a list that will host the results of each model
for(i in 1:length(levels(dat$Study_ID))){ # N models = N studies
VCV_matrix[[i]]<-make_VCV_matrix(dat[dat$Study_ID != levels(dat$Study_ID)[i], ], V="var_lnRR", cluster="Cohort_ID", obs="Effect_ID") # Create a new VCV matrix for each new model
Leave1studyout[[i]] <- rma.mv(yi = lnRR, V = VCV_matrix[[i]], # Same model structure as all the models we fitted
random = list(~1|Study_ID,
~1|Phylogeny,
~1|Species_common,
~1|PFAS_type,
~1|Effect_ID),
R= list(Phylogeny = cor_tree),
test = "t",
data = dat[dat$Study_ID != levels(dat$Study_ID)[i], ]) # Generate a new model for each new data (iterative removal of one study at a time)
}
# The output is a list so we need to summarise the coefficients of all the models performed
results.Leave1studyout<-as.data.frame(cbind(
sapply(Leave1studyout, function(x) summary(x)$beta), # extract the beta coefficient from all models
sapply(Leave1studyout, function(x) summary(x)$se), # extract the standard error from all models
sapply(Leave1studyout, function(x) summary(x)$zval), # extract the z value from all models
sapply(Leave1studyout, function(x) summary(x)$pval), # extract the p value from all models
sapply(Leave1studyout, function(x) summary(x)$ci.lb), # extract the lower confidence interval for all models
sapply(Leave1studyout, function(x) summary(x)$ci.ub))) # extract the upper confidence interval for all models
colnames(results.Leave1studyout)=c("Estimate", "SE", "zval", "pval", "ci.lb", "ci.ub") # change column names
kable(results.Leave1studyout)%>% kable_styling("striped", position="left") %>% scroll_box(width="100%", height="500px") # Table of the results from all models| Estimate | SE | zval | pval | ci.lb | ci.ub |
|---|---|---|---|---|---|
| -0.3253221 | 0.3107507 | -1.0468911 | 0.2956467 | -0.9358339 | 0.2851897 |
| -0.4037084 | 0.3101145 | -1.3018043 | 0.1935803 | -1.0129906 | 0.2055738 |
| -0.3997524 | 0.3468279 | -1.1525957 | 0.2497971 | -1.0816832 | 0.2821784 |
| 0.0435382 | 0.2731984 | 0.1593648 | 0.8734478 | -0.4932604 | 0.5803368 |
| -0.3312637 | 0.3129344 | -1.0585724 | 0.2903281 | -0.9461575 | 0.2836301 |
| -0.2423434 | 0.3010346 | -0.8050351 | 0.4211923 | -0.8338304 | 0.3491436 |
| -0.3309747 | 0.3124412 | -1.0593181 | 0.2899684 | -0.9448401 | 0.2828908 |
| -0.2229376 | 0.3086359 | -0.7223322 | 0.4706194 | -0.8301566 | 0.3842813 |
| -0.3882687 | 0.3207704 | -1.2104253 | 0.2267090 | -1.0185498 | 0.2420125 |
| -0.4843182 | 0.2868896 | -1.6881692 | 0.0920640 | -1.0481112 | 0.0794748 |
## # A tibble: 10 x 3
## # Groups: Author_year [10]
## Author_year Study_ID mean
## <chr> <fct> <dbl>
## 1 Alves_2017 F001 -0.0774
## 2 Barbosa_2018 F002 0.198
## 3 Bhavsar_2014 F003 0.153
## 4 DelGobbo_2008 F005 -2.00
## 5 Hu_2020 F006 -0.134
## 6 Kim_2020 F007 -0.887
## 7 Luo_2019 F008 -0.161
## 8 Sungur_2019 F010 -0.893
## 9 Taylor_2019 F011 0.213
## 10 Vassiliadou_2015 F013 0.673
Study_ID F005 (Del Gobbo et al. 2008)dat.sens<-filter(dat, Author_year!="DelGobbo_2008")
include <- row.names(cor_tree) %in% dat.sens$Phylogeny # Check which rows are present in the phylogenetic tree
cor_tree_sens <- cor_tree[include, include] # Only include the species that match the reduced data set
dat.sens<-as.data.frame(dat.sens) # convert data set into a data frame to calculate VCV matrix
VCV_lnRR.sens<-make_VCV_matrix(dat.sens, V = "var_lnRR", cluster = "Cohort_ID", obs = "Effect_ID", rho = 0.5) # create VCV matrix for the specified data
mod.sens<- rma.mv(lnRR, VCV_lnRR.sens,
mods=~Length_cooking_time_in_s,
random = list(~1|Study_ID,
~1|Phylogeny,
~1|Species_common,
~1|PFAS_type,
~1|Effect_ID),
R= list(Phylogeny = cor_tree_sens),
test = "t",
data = dat.sens)
summary(mod.sens)##
## Multivariate Meta-Analysis Model (k = 430; method: REML)
##
## logLik Deviance AIC BIC AICc
## -276.6805 553.3611 567.3611 595.7749 567.6277
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.2023 0.4498 8 no Study_ID no
## sigma^2.2 0.0311 0.1763 22 no Phylogeny yes
## sigma^2.3 0.0112 0.1059 22 no Species_common no
## sigma^2.4 0.0964 0.3105 17 no PFAS_type no
## sigma^2.5 0.0683 0.2613 430 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 428) = 1249.4809, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 428) = 83.0392, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.6929 0.2349 2.9499 0.0034 0.2312 1.1546
## Length_cooking_time_in_s -0.0012 0.0001 -9.1126 <.0001 -0.0015 -0.0009
##
## intrcpt **
## Length_cooking_time_in_s ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dat.time.sens<-filter(dat.sens, Length_cooking_time_in_s!="NA")
plot_continuous(dat.time.sens, mod.sens, dat.time.sens$Length_cooking_time_in_s, "Cooking time (s)") # The relationship with cooking time appears even stronger oil_dat.sens <- filter(dat.sens, Cooking_Category=="oil-based")
water_dat.sens <- filter(dat.sens, Cooking_Category=="water-based")
dry_dat.sens <- filter(dat.sens, Cooking_Category=="No liquid")
oil_dat_time.sens<-filter(oil_dat.sens, Length_cooking_time_in_s!="NA")
water_dat_time.sens<-filter(water_dat.sens, Length_cooking_time_in_s!="NA")
dry_dat_time.sens<-filter(dry_dat.sens, Length_cooking_time_in_s!="NA")
model_oil_time.sens<-run_model_oil(oil_dat_time.sens, ~Length_cooking_time_in_s)
model_water_time.sens<-run_model_water(water_dat_time.sens, ~Length_cooking_time_in_s)
model_dry_time.sens<-run_model_dry(dry_dat_time.sens, ~Length_cooking_time_in_s)
pred_oil_time.sens<-predict.rma(model_oil_time.sens)
pred_water_time.sens<-predict.rma(model_water_time.sens)
pred_dry_time.sens<-predict.rma(model_dry_time.sens)
oil_dat_time.sens<-mutate(oil_dat_time.sens,
ci.lb = pred_oil_time.sens$ci.lb,
ci.ub = pred_oil_time.sens$ci.ub,
fit = pred_oil_time.sens$pred)
water_dat_time.sens<-mutate(water_dat_time.sens,
ci.lb = pred_water_time.sens$ci.lb,
ci.ub = pred_water_time.sens$ci.ub,
fit = pred_water_time.sens$pred)
dry_dat_time.sens<-mutate(dry_dat_time.sens,
ci.lb = pred_dry_time.sens$ci.lb,
ci.ub = pred_dry_time.sens$ci.ub,
fit = pred_dry_time.sens$pred)
# Actual plot
ggplot(dat.sens,aes(x = Length_cooking_time_in_s, y = lnRR, fill=Cooking_Category)) +
geom_ribbon(data=water_dat_time.sens, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=water_dat_time.sens,aes(y = fit), size = 1.5, col="dodgerblue")+
geom_ribbon(data=oil_dat_time.sens, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.3) +
geom_line(data=oil_dat_time.sens,aes(y = fit), size = 1.5, col="goldenrod")+
geom_ribbon(data=dry_dat_time.sens, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.25) +
geom_line(data=dry_dat_time.sens,aes(y = fit), size = 1.5, col="palegreen3")+
geom_point(aes(size=(1/sqrt(var_lnRR)), fill=Cooking_Category), shape=21, alpha=0.8) +
scale_fill_manual(values=c("#55C667FF", "goldenrod2", "dodgerblue3"))+
labs(x = "Cooking time (s)", y = "lnRR", size = "Precison (1/SE)") +
scale_size_continuous(range=c(1,10))+
theme_bw() +
geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+
theme(text = element_text(size = 18, colour = "black", hjust = 0.5),
legend.text=element_text(size=14),
legend.position=c(0,0),
legend.justification = c(0,0),
legend.background = element_blank(),
legend.direction="horizontal",
legend.title = element_text(size=15),
panel.border=element_rect(colour="black", fill=NA, size=1.2))dat.sens.vol<-filter(dat.sens, Volume_liquid_ml!="NA")
include <- row.names(cor_tree) %in% dat.sens.vol$Phylogeny # Check which rows are present in the phylogenetic tree
cor_tree_sens.vol <- cor_tree[include, include] # Only include the species that match the reduced data set
VCV_lnRR.sens.vol<-make_VCV_matrix(dat.sens.vol, V = "var_lnRR", cluster = "Cohort_ID", obs = "Effect_ID", rho = 0.5) # create VCV matrix for the specified data
mod.sens.vol<- rma.mv(lnRR, VCV_lnRR.sens.vol,
mods=~ log(Volume_liquid_ml),
random = list(~1|Study_ID,
~1|Phylogeny,
~1|Species_common,
~1|PFAS_type,
~1|Effect_ID),
R= list(Phylogeny = cor_tree_sens.vol),
test = "t",
data = dat.sens.vol)
summary(mod.sens.vol)##
## Multivariate Meta-Analysis Model (k = 413; method: REML)
##
## logLik Deviance AIC BIC AICc
## -389.4527 778.9053 792.9053 821.0355 793.1832
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.3293 0.5739 7 no Study_ID no
## sigma^2.2 0.0413 0.2031 26 no Phylogeny yes
## sigma^2.3 0.1058 0.3252 27 no Species_common no
## sigma^2.4 0.1297 0.3601 18 no PFAS_type no
## sigma^2.5 0.2120 0.4604 413 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 411) = 3390.2773, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 411) = 0.4459, p-val = 0.5046
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -0.3392 0.4516 -0.7511 0.4530 -1.2268 0.5485
## log(Volume_liquid_ml) 0.0462 0.0691 0.6678 0.5046 -0.0897 0.1821
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The effect of volume of liquid is entirely driven by one study!!
oil_dat.sens <- filter(dat.sens, Cooking_Category=="oil-based")
water_dat.sens <- filter(dat.sens, Cooking_Category=="water-based")
dry_dat.sens <- filter(dat.sens, Cooking_Category=="No liquid")
oil_dat_vol.sens<-filter(oil_dat.sens, Volume_liquid_ml!="NA")
water_dat_vol.sens<-filter(water_dat.sens, Volume_liquid_ml!="NA")
dry_dat_vol.sens<-filter(dry_dat.sens, Volume_liquid_ml!="NA")
model_oil_vol.sens<-run_model_oil(oil_dat_vol.sens, ~log(Volume_liquid_ml))
model_water_vol.sens<-run_model_water(water_dat_vol.sens, ~log(Volume_liquid_ml))
model_dry_vol.sens<-run_model_dry(dry_dat_vol.sens, ~log(Volume_liquid_ml))
pred_oil_vol.sens<-predict.rma(model_oil_vol.sens)
pred_water_vol.sens<-predict.rma(model_water_vol.sens)
pred_dry_vol.sens<-predict.rma(model_dry_vol.sens)
oil_dat_vol.sens<-mutate(oil_dat_vol.sens,
ci.lb = pred_oil_vol.sens$ci.lb,
ci.ub = pred_oil_vol.sens$ci.ub,
fit = pred_oil_vol.sens$pred)
water_dat_vol.sens<-mutate(water_dat_vol.sens,
ci.lb = pred_water_vol.sens$ci.lb,
ci.ub = pred_water_vol.sens$ci.ub,
fit = pred_water_vol.sens$pred)
dry_dat_vol.sens<-mutate(dry_dat_vol.sens,
ci.lb = pred_dry_vol.sens$ci.lb,
ci.ub = pred_dry_vol.sens$ci.ub,
fit = pred_dry_vol.sens$pred)
# Actual plot
ggplot(dat.sens,aes(x = log(Volume_liquid_ml), y = lnRR, fill=Cooking_Category)) +
geom_ribbon(data=water_dat_vol.sens, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=water_dat_vol.sens,aes(y = fit), size = 1.5, col="dodgerblue")+
geom_ribbon(data=oil_dat_vol.sens, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.3) +
geom_line(data=oil_dat_vol.sens,aes(y = fit), size = 1.5, col="goldenrod")+
geom_ribbon(data=dry_dat_vol.sens, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=dry_dat_vol.sens,aes(y = fit), size = 1.5, col="palegreen3")+
geom_point(aes(size=(1/sqrt(var_lnRR)), fill=Cooking_Category), shape=21, alpha=0.8) +
scale_fill_manual(values=c("#55C667FF", "goldenrod2", "dodgerblue3"))+
labs(x = "ln(Volume of liquid (mL))", y = "lnRR", size = "Precison (1/SE)") +
scale_size_continuous(range=c(1,10))+
theme_bw() +
geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+
theme(text = element_text(size = 18, colour = "black", hjust = 0.5),
legend.text=element_text(size=14),
legend.position=c(0,0),
legend.justification = c(0,0),
legend.background = element_blank(),
legend.direction="horizontal",
legend.title = element_text(size=15),
panel.border=element_rect(colour="black", fill=NA, size=1.2))dat.sens.PFAS<-filter(dat.sens, PFAS_carbon_chain!="NA")
include <- row.names(cor_tree) %in% dat.sens.PFAS$Phylogeny # Check which rows are present in the phylogenetic tree
cor_tree_sens.PFAS <- cor_tree[include, include] # Only include the species that match the reduced data set
VCV_lnRR.sens.PFAS<-make_VCV_matrix(dat.sens.PFAS, V = "var_lnRR", cluster = "Cohort_ID", obs = "Effect_ID", rho = 0.5) # create VCV matrix for the specified data
mod.sens.PFAS<- rma.mv(lnRR, VCV_lnRR.sens.PFAS,
mods=~ PFAS_carbon_chain,
random = list(~1|Study_ID,
~1|Phylogeny,
~1|Species_common,
~1|PFAS_type,
~1|Effect_ID),
R= list(Phylogeny = cor_tree_sens.PFAS),
test = "t",
data = dat.sens.PFAS)
summary(mod.sens.PFAS)##
## Multivariate Meta-Analysis Model (k = 486; method: REML)
##
## logLik Deviance AIC BIC AICc
## -464.4535 928.9070 942.9070 972.1816 943.1423
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.2361 0.4859 9 no Study_ID no
## sigma^2.2 0.0998 0.3159 30 no Phylogeny yes
## sigma^2.3 0.1036 0.3218 31 no Species_common no
## sigma^2.4 0.0968 0.3111 18 no PFAS_type no
## sigma^2.5 0.2299 0.4795 486 no Effect_ID no
##
## Test for Residual Heterogeneity:
## QE(df = 484) = 4439.7079, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 484) = 1.1682, p-val = 0.2803
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -0.2244 0.3706 -0.6054 0.5452 -0.9525 0.5038
## PFAS_carbon_chain 0.0301 0.0278 1.0809 0.2803 -0.0246 0.0847
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
oil_dat.sens <- filter(dat.sens, Cooking_Category=="oil-based")
water_dat.sens <- filter(dat.sens, Cooking_Category=="water-based")
dry_dat.sens <- filter(dat.sens, Cooking_Category=="No liquid")
oil_dat_PFAS.sens<-filter(oil_dat.sens, PFAS_carbon_chain!="NA")
water_dat_PFAS.sens<-filter(water_dat.sens, PFAS_carbon_chain!="NA")
dry_dat_PFAS.sens<-filter(dry_dat.sens, PFAS_carbon_chain!="NA")
model_oil_PFAS.sens<-run_model_oil(oil_dat_PFAS.sens, ~PFAS_carbon_chain)
model_water_PFAS.sens<-run_model_water(water_dat_PFAS.sens, ~PFAS_carbon_chain)
model_dry_PFAS.sens<-run_model_dry(dry_dat_PFAS.sens, ~PFAS_carbon_chain)
pred_oil_PFAS.sens<-predict.rma(model_oil_PFAS.sens)
pred_water_PFAS.sens<-predict.rma(model_water_PFAS.sens)
pred_dry_PFAS.sens<-predict.rma(model_dry_PFAS.sens)
oil_dat_PFAS.sens<-mutate(oil_dat_PFAS.sens,
ci.lb = pred_oil_PFAS.sens$ci.lb,
ci.ub = pred_oil_PFAS.sens$ci.ub,
fit = pred_oil_PFAS.sens$pred)
water_dat_PFAS.sens<-mutate(water_dat_PFAS.sens,
ci.lb = pred_water_PFAS.sens$ci.lb,
ci.ub = pred_water_PFAS.sens$ci.ub,
fit = pred_water_PFAS.sens$pred)
dry_dat_PFAS.sens<-mutate(dry_dat_PFAS.sens,
ci.lb = pred_dry_PFAS.sens$ci.lb,
ci.ub = pred_dry_PFAS.sens$ci.ub,
fit = pred_dry_PFAS.sens$pred)
# Actual plot
ggplot(dat.sens,aes(x = PFAS_carbon_chain, y = lnRR, fill=Cooking_Category)) +
geom_ribbon(data=dry_dat_PFAS.sens, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.2) +
geom_line(data=dry_dat_PFAS.sens,aes(y = fit), size = 1.5, col="palegreen3")+
geom_ribbon(data=oil_dat_PFAS.sens, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.3) +
geom_line(data=oil_dat_PFAS.sens,aes(y = fit), size = 1.5, col="goldenrod")+
geom_ribbon(data=water_dat_PFAS.sens, aes(ymin = ci.lb, ymax = ci.ub, color = NULL), alpha = 0.3) +
geom_line(data=water_dat_PFAS.sens,aes(y = fit), size = 1.5, col="dodgerblue")+
geom_point(aes(size=(1/sqrt(var_lnRR)), fill=Cooking_Category), shape=21, alpha=0.8) +
scale_fill_manual(values=c("#55C667FF", "goldenrod2", "dodgerblue3"))+
labs(x = "PFAS carbon chain length", y = "lnRR", size = "Precison (1/SE)") +
scale_size_continuous(range=c(1,10))+
theme_bw() +
geom_hline(yintercept = 0,linetype = 2, colour = "black",alpha=0.5)+
theme(text = element_text(size = 18, colour = "black", hjust = 0.5),
legend.text=element_text(size=14),
legend.position=c(0,0),
legend.justification = c(0,0),
legend.background = element_blank(),
legend.direction="horizontal",
legend.title = element_text(size=15),
panel.border=element_rect(colour="black", fill=NA, size=1.2))